

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>agpy.psf_fitter &mdash; agpy 0.1.2 documentation</title>
    
    <link rel="stylesheet" href="../../_static/extra.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '0.1.2',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <link rel="top" title="agpy 0.1.2 documentation" href="../../index.html" />
    <link rel="up" title="agpy" href="../agpy.html" />
     
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setDomainName', 'pyspeckit.bitbucket.org']);
      _gaq.push(['_setAllowHash', false]);
      _gaq.push(['_trackPageview']);


    </script>
    <link rel="stylesheet" type="text/css" href="../../_static/extra.css" />
  </head>
  <body>
    <div class="header-wrapper">
      <div class="header">
        <h1><a href="../../index.html">agpy 0.1.2 documentation</a></h1>
        <div class="rel">
          <a href="http://agpy.googlecode.com">agpy Home </a>  |
          <a href=../../index.html>Docs Home </a>  |
          <a href="http://code.google.com/p/agpy/w/list">Wiki</a>  |
          <a href=../../search.html>Search </a>
        </div>
       </div>
    </div>

    <div class="content-wrapper">
      <div class="content">
        <div class="document">
            
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <h1>Source code for agpy.psf_fitter</h1><div class="highlight"><pre>
<span class="c"># Fit a PSF of type Airy or Gaussian...</span>
<span class="kn">from</span> <span class="nn">gaussfitter</span> <span class="kn">import</span> <span class="n">twodgaussian</span><span class="p">,</span><span class="n">moments</span>
<span class="kn">import</span> <span class="nn">numpy</span>
<span class="kn">import</span> <span class="nn">scipy</span>
<span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="n">pi</span>
<span class="kn">from</span> <span class="nn">agpy.mpfit</span> <span class="kn">import</span> <span class="n">mpfit</span>

<span class="k">def</span> <span class="nf">_airy_func</span><span class="p">(</span><span class="n">rr</span><span class="p">,</span> <span class="n">amplitude</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">width</span><span class="o">=</span><span class="mf">1.0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    For a simple radially symmetric airy function, returns the value at a given</span>
<span class="sd">    (normalized) radius</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="n">amplitude</span> <span class="o">*</span> <span class="p">(</span><span class="mf">2.0</span> <span class="o">*</span> <span class="n">scipy</span><span class="o">.</span><span class="n">special</span><span class="o">.</span><span class="n">j1</span><span class="p">(</span><span class="n">rr</span><span class="o">/</span><span class="n">width</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">rr</span><span class="o">/</span><span class="n">width</span><span class="p">))</span><span class="o">**</span><span class="mi">2</span>

<span class="k">def</span> <span class="nf">_gaussian_func</span><span class="p">(</span><span class="n">rr</span><span class="p">,</span> <span class="n">amplitude</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">sigma</span><span class="o">=</span><span class="mf">1.0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    For a simple radially symmetric Gaussian function, returns the value at a given</span>
<span class="sd">    (normalized) radius</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="n">amplitude</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">rr</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="mf">2.0</span> <span class="o">*</span> <span class="n">sigma</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="p">)</span>

<div class="viewcode-block" id="airy"><a class="viewcode-back" href="../../agpy.html#agpy.psf_fitter.airy">[docs]</a><span class="k">def</span> <span class="nf">airy</span><span class="p">(</span><span class="n">inpars</span><span class="p">,</span> <span class="n">circle</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span> <span class="n">rotate</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span> <span class="n">vheight</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span> <span class="n">shape</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">fwhm</span><span class="o">=</span><span class="bp">False</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Returns a 2d Airy *function* of the form:</span>
<span class="sd">        x&#39; = numpy.cos(rota) * x - numpy.sin(rota) * y</span>
<span class="sd">        y&#39; = numpy.sin(rota) * x + numpy.cos(rota) * y</span>
<span class="sd">        (rota should be in degrees)</span>
<span class="sd">        radius = sqrt( (x&#39;-xcen)^2 + (y&#39;-ycen)^2 )</span>
<span class="sd">        g = b + a * 2.0*BesselJ1( radius ) / radius</span>
<span class="sd">            (with a correction for the divide-by-zero in the center)</span>

<span class="sd">        inpars = [b,a,center_x,center_y,width_x,width_y,rota]</span>
<span class="sd">                 (b is background height, a is peak amplitude)</span>

<span class="sd">        where x and y are the input parameters of the returned function,</span>
<span class="sd">        and all other parameters are specified by this function</span>

<span class="sd">        However, the above values are passed by list.  The list should be:</span>
<span class="sd">        inpars = (height,amplitude,center_x,center_y,width_x,width_y,rota)</span>

<span class="sd">        You can choose to ignore / neglect some of the above input parameters </span>
<span class="sd">            unumpy.sing the following options:</span>
<span class="sd">            circle=1 - default is a circular Airy Disk.  An elliptical Airy is</span>
<span class="sd">                possible, but probably not physically motivated (unless it&#39;s</span>
<span class="sd">                sampled onto a stretched grid?).</span>
<span class="sd">            rotate=0 - default allows rotation of the gaussian ellipse.  Can</span>
<span class="sd">                remove last parameter by setting rotate=0</span>
<span class="sd">            vheight=1 - default allows a variable height-above-zero, i.e. an</span>
<span class="sd">                additive constant for the Airy function.  Can remove first</span>
<span class="sd">                parameter by setting this to 0</span>
<span class="sd">            shape=None - if shape is set (to a 2-parameter list) then returns</span>
<span class="sd">                an image with the gaussian defined by inpars</span>
<span class="sd">        fwhm - if set, assumes the Width parameters input are FWHM widths, so</span>
<span class="sd">            they&#39;ll be converted to &quot;Sigma&quot; widths by s = FWHM/2.0/1.61633</span>
<span class="sd">            (http://en.wikipedia.org/wiki/Airy_disk </span>
<span class="sd">            and http://home.fnal.gov/~neilsen/notebook/astroPSF/astroPSF.html)</span>
<span class="sd">        &quot;&quot;&quot;</span>
    <span class="n">inpars_old</span> <span class="o">=</span> <span class="n">inpars</span>
    <span class="n">inpars</span> <span class="o">=</span> <span class="nb">list</span><span class="p">(</span><span class="n">inpars</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">vheight</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">height</span> <span class="o">=</span> <span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
        <span class="n">height</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">height</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">height</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">amplitude</span><span class="p">,</span> <span class="n">center_y</span><span class="p">,</span> <span class="n">center_x</span> <span class="o">=</span> <span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">),</span><span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">),</span><span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">amplitude</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">amplitude</span><span class="p">)</span>
    <span class="n">center_x</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">center_x</span><span class="p">)</span>
    <span class="n">center_y</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">center_y</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">circle</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">width</span> <span class="o">=</span> <span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
        <span class="n">width_x</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">width</span><span class="p">)</span>
        <span class="n">width_y</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">width</span><span class="p">)</span>
        <span class="n">rotate</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">width_x</span><span class="p">,</span> <span class="n">width_y</span> <span class="o">=</span> <span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">),</span><span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
        <span class="n">width_x</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">width_x</span><span class="p">)</span>
        <span class="n">width_y</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">width_y</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">rotate</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">rota</span> <span class="o">=</span> <span class="n">inpars</span><span class="o">.</span><span class="n">pop</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
        <span class="n">rota</span> <span class="o">=</span> <span class="n">pi</span><span class="o">/</span><span class="mf">180.</span> <span class="o">*</span> <span class="nb">float</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span>
        <span class="n">rcen_x</span> <span class="o">=</span> <span class="n">center_x</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span> <span class="o">-</span> <span class="n">center_y</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span>
        <span class="n">rcen_y</span> <span class="o">=</span> <span class="n">center_x</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span> <span class="o">+</span> <span class="n">center_y</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">rcen_x</span> <span class="o">=</span> <span class="n">center_x</span>
        <span class="n">rcen_y</span> <span class="o">=</span> <span class="n">center_y</span>
    <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">inpars</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&quot;There are still input parameters:&quot;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">inpars</span><span class="p">)</span> <span class="o">+</span> \
                <span class="s">&quot; and you&#39;ve input: &quot;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">inpars_old</span><span class="p">)</span> <span class="o">+</span> \
                <span class="s">&quot; circle=</span><span class="si">%d</span><span class="s">, rotate=</span><span class="si">%d</span><span class="s">, vheight=</span><span class="si">%d</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">)</span> <span class="p">)</span>

    <span class="k">if</span> <span class="n">fwhm</span><span class="p">:</span>
        <span class="n">width_x</span> <span class="o">/=</span> <span class="mf">2.0</span> <span class="o">*</span> <span class="mf">1.61633</span>
        <span class="n">width_y</span> <span class="o">/=</span> <span class="mf">2.0</span> <span class="o">*</span> <span class="mf">1.61633</span>
            
    <span class="k">def</span> <span class="nf">rotairy</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">rotate</span><span class="o">==</span><span class="mi">1</span><span class="p">:</span>
            <span class="n">xp</span> <span class="o">=</span> <span class="n">x</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span> <span class="o">-</span> <span class="n">y</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span>
            <span class="n">yp</span> <span class="o">=</span> <span class="n">x</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span> <span class="o">+</span> <span class="n">y</span> <span class="o">*</span> <span class="n">numpy</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">rota</span><span class="p">)</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">xp</span> <span class="o">=</span> <span class="n">x</span>
            <span class="n">yp</span> <span class="o">=</span> <span class="n">y</span>
        <span class="n">rr</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(((</span><span class="n">rcen_x</span><span class="o">-</span><span class="n">xp</span><span class="p">)</span><span class="o">/</span><span class="n">width_x</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span>
             <span class="p">((</span><span class="n">rcen_y</span><span class="o">-</span><span class="n">yp</span><span class="p">)</span><span class="o">/</span><span class="n">width_y</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        <span class="c"># http://en.wikipedia.org/wiki/Airy_disk</span>
        <span class="n">airy_func</span> <span class="o">=</span> <span class="p">(</span><span class="mf">2.0</span> <span class="o">*</span> <span class="n">scipy</span><span class="o">.</span><span class="n">special</span><span class="o">.</span><span class="n">j1</span><span class="p">(</span><span class="n">rr</span><span class="p">)</span> <span class="o">/</span> <span class="n">rr</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span>
        <span class="n">airy_func</span><span class="p">[</span><span class="n">rr</span><span class="o">==</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
        <span class="n">airy</span> <span class="o">=</span> <span class="n">height</span> <span class="o">+</span> <span class="n">amplitude</span> <span class="o">*</span> <span class="n">airy_func</span>

        <span class="k">return</span> <span class="n">airy</span>
    <span class="k">if</span> <span class="n">shape</span> <span class="ow">is</span> <span class="ow">not</span> <span class="bp">None</span><span class="p">:</span>
        <span class="k">return</span> <span class="n">rotairy</span><span class="p">(</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">indices</span><span class="p">(</span><span class="n">shape</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="k">return</span> <span class="n">rotairy</span>
</div>
<div class="viewcode-block" id="psffit"><a class="viewcode-back" href="../../agpy.html#agpy.psf_fitter.psffit">[docs]</a><span class="k">def</span> <span class="nf">psffit</span><span class="p">(</span><span class="n">data</span><span class="p">,</span><span class="n">err</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span><span class="n">params</span><span class="o">=</span><span class="p">[],</span><span class="n">autoderiv</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span><span class="n">return_all</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span><span class="n">circle</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span>
        <span class="n">fixed</span><span class="o">=</span><span class="n">numpy</span><span class="o">.</span><span class="n">repeat</span><span class="p">(</span><span class="bp">False</span><span class="p">,</span><span class="mi">7</span><span class="p">),</span><span class="n">limitedmin</span><span class="o">=</span><span class="p">[</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">True</span><span class="p">,</span><span class="bp">True</span><span class="p">,</span><span class="bp">True</span><span class="p">],</span>
        <span class="n">limitedmax</span><span class="o">=</span><span class="p">[</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">False</span><span class="p">,</span><span class="bp">True</span><span class="p">],</span>
        <span class="n">usemoment</span><span class="o">=</span><span class="n">numpy</span><span class="o">.</span><span class="n">array</span><span class="p">([],</span><span class="n">dtype</span><span class="o">=</span><span class="s">&#39;bool&#39;</span><span class="p">),</span>
        <span class="n">minpars</span><span class="o">=</span><span class="n">numpy</span><span class="o">.</span><span class="n">repeat</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">7</span><span class="p">),</span><span class="n">maxpars</span><span class="o">=</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">360</span><span class="p">],</span>
        <span class="n">rotate</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">vheight</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span><span class="n">quiet</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span><span class="n">returnmp</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span>
        <span class="n">returnfitimage</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span>
        <span class="n">psffunction</span><span class="o">=</span><span class="n">airy</span><span class="p">,</span> 
        <span class="n">extra_pars</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span>
        <span class="n">return_parinfo</span><span class="o">=</span><span class="bp">False</span><span class="p">,</span>
        <span class="o">**</span><span class="n">kwargs</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    PSF fitter with the ability to fit a variety of different forms of</span>
<span class="sd">    2-dimensional gaussian OR an Airy.</span>
<span class="sd">    This code is mostly directly copied from gaussfitter.py and presents </span>
<span class="sd">    yet another argument for me turning this into a class...</span>
<span class="sd">    </span>
<span class="sd">    Input Parameters:</span>
<span class="sd">        data - 2-dimensional data array</span>
<span class="sd">        err=None - error array with same size as data array</span>
<span class="sd">        params=[] - initial input parameters for Gaussian function.</span>
<span class="sd">            (height, amplitude, x, y, width_x, width_y, rota)</span>
<span class="sd">            if not input, these will be determined from the moments of the system, </span>
<span class="sd">            assuming no rotation</span>
<span class="sd">        autoderiv=1 - use the autoderiv provided in the lmder.f function (the</span>
<span class="sd">            alternative is to us an analytic derivative with lmdif.f: this method</span>
<span class="sd">            is less robust)</span>
<span class="sd">        return_all=0 - Default is to return only the Gaussian parameters.  </span>
<span class="sd">                   1 - fit params, fit error</span>
<span class="sd">        returnfitimage - returns (best fit params,best fit image)</span>
<span class="sd">        returnmp - returns the full mpfit struct</span>
<span class="sd">        circle=0 - default is an elliptical gaussian (different x, y widths),</span>
<span class="sd">            but can reduce the input by one parameter if it&#39;s a circular gaussian</span>
<span class="sd">        rotate=1 - default allows rotation of the gaussian ellipse.  Can remove</span>
<span class="sd">            last parameter by setting rotate=0.  numpy.expects angle in DEGREES</span>
<span class="sd">        vheight=1 - default allows a variable height-above-zero, i.e. an</span>
<span class="sd">            additive constant for the Gaussian function.  Can remove first</span>
<span class="sd">            parameter by setting this to 0</span>
<span class="sd">        usemoment - can choose which parameters to use a moment estimation for.</span>
<span class="sd">            Other parameters will be taken from params.  Needs to be a boolean</span>
<span class="sd">            array.</span>
<span class="sd">        extra_pars - If your psffunction requires extra parameters, pass their</span>
<span class="sd">            parinfo dictionaries through this variable</span>

<span class="sd">    Output:</span>
<span class="sd">        Default output is a set of Gaussian parameters with the same shape as</span>
<span class="sd">            the input parameters</span>

<span class="sd">        Can also output the covariance matrix, &#39;infodict&#39; that contains a lot</span>
<span class="sd">            more detail about the fit (see scipy.optimize.leastsq), and a message</span>
<span class="sd">            from leastsq telling what the exit status of the fitting routine was</span>

<span class="sd">        Warning: Does NOT necessarily output a rotation angle between 0 and 360 degrees.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">usemoment</span><span class="o">=</span><span class="n">numpy</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">usemoment</span><span class="p">,</span><span class="n">dtype</span><span class="o">=</span><span class="s">&#39;bool&#39;</span><span class="p">)</span>
    <span class="n">params</span><span class="o">=</span><span class="n">numpy</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">params</span><span class="p">,</span><span class="n">dtype</span><span class="o">=</span><span class="s">&#39;float&#39;</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">usemoment</span><span class="o">.</span><span class="n">any</span><span class="p">()</span> <span class="ow">and</span> <span class="nb">len</span><span class="p">(</span><span class="n">params</span><span class="p">)</span><span class="o">==</span><span class="nb">len</span><span class="p">(</span><span class="n">usemoment</span><span class="p">):</span>
        <span class="n">moment</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">moments</span><span class="p">(</span><span class="n">data</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">,</span><span class="o">**</span><span class="n">kwargs</span><span class="p">),</span><span class="n">dtype</span><span class="o">=</span><span class="s">&#39;float&#39;</span><span class="p">)</span>
        <span class="n">params</span><span class="p">[</span><span class="n">usemoment</span><span class="p">]</span> <span class="o">=</span> <span class="n">moment</span><span class="p">[</span><span class="n">usemoment</span><span class="p">]</span>
    <span class="k">elif</span> <span class="n">params</span> <span class="o">==</span> <span class="p">[]</span> <span class="ow">or</span> <span class="nb">len</span><span class="p">(</span><span class="n">params</span><span class="p">)</span><span class="o">==</span><span class="mi">0</span><span class="p">:</span>
        <span class="n">params</span> <span class="o">=</span> <span class="p">(</span><span class="n">moments</span><span class="p">(</span><span class="n">data</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">,</span><span class="o">**</span><span class="n">kwargs</span><span class="p">))</span>
    <span class="k">if</span> <span class="n">vheight</span><span class="o">==</span><span class="mi">0</span><span class="p">:</span>
        <span class="n">vheight</span><span class="o">=</span><span class="mi">1</span>
        <span class="n">params</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">concatenate</span><span class="p">([[</span><span class="mi">0</span><span class="p">],</span><span class="n">params</span><span class="p">])</span>
        <span class="n">fixed</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mi">1</span>


    <span class="c"># mpfit will fail if it is given a start parameter outside the allowed range:</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">params</span><span class="p">)):</span> 
        <span class="k">if</span> <span class="n">params</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">&gt;</span> <span class="n">maxpars</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="ow">and</span> <span class="n">limitedmax</span><span class="p">[</span><span class="n">i</span><span class="p">]:</span> <span class="n">params</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">maxpars</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
        <span class="k">if</span> <span class="n">params</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">minpars</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="ow">and</span> <span class="n">limitedmin</span><span class="p">[</span><span class="n">i</span><span class="p">]:</span> <span class="n">params</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">minpars</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>

    <span class="k">if</span> <span class="n">err</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
        <span class="n">errorfunction</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">p</span><span class="p">:</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ravel</span><span class="p">((</span><span class="n">psffunction</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">)</span>\
                <span class="p">(</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">indices</span><span class="p">(</span><span class="n">data</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span> <span class="o">-</span> <span class="n">data</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">errorfunction</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">p</span><span class="p">:</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ravel</span><span class="p">((</span><span class="n">psffunction</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">)</span>\
                <span class="p">(</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">indices</span><span class="p">(</span><span class="n">data</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span> <span class="o">-</span> <span class="n">data</span><span class="p">)</span><span class="o">/</span><span class="n">err</span><span class="p">)</span>
    <span class="k">def</span> <span class="nf">mpfitfun</span><span class="p">(</span><span class="n">data</span><span class="p">,</span><span class="n">err</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">err</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
            <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">fjac</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span> <span class="k">return</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">numpy</span><span class="o">.</span><span class="n">ravel</span><span class="p">(</span><span class="n">data</span><span class="o">-</span><span class="n">psffunction</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">)</span>\
                    <span class="p">(</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">indices</span><span class="p">(</span><span class="n">data</span><span class="o">.</span><span class="n">shape</span><span class="p">)))]</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">fjac</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span> <span class="k">return</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="n">numpy</span><span class="o">.</span><span class="n">ravel</span><span class="p">((</span><span class="n">data</span><span class="o">-</span><span class="n">psffunction</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">)</span>\
                    <span class="p">(</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">indices</span><span class="p">(</span><span class="n">data</span><span class="o">.</span><span class="n">shape</span><span class="p">)))</span><span class="o">/</span><span class="n">err</span><span class="p">)]</span>
        <span class="k">return</span> <span class="n">f</span>

                    
    <span class="n">parinfo</span> <span class="o">=</span> <span class="p">[</span> 
                <span class="p">{</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">1</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">1</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">1</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;AMPLITUDE&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">},</span>
                <span class="p">{</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">2</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">2</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">2</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;XSHIFT&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">},</span>
                <span class="p">{</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">3</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">3</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">3</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;YSHIFT&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">},</span>
                <span class="p">{</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">4</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">4</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">4</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;XWIDTH&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">}</span> <span class="p">]</span>
    <span class="k">if</span> <span class="n">vheight</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="n">parinfo</span><span class="o">.</span><span class="n">insert</span><span class="p">(</span><span class="mi">0</span><span class="p">,{</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">0</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">0</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;HEIGHT&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">})</span>
    <span class="k">if</span> <span class="n">circle</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">parinfo</span><span class="o">.</span><span class="n">append</span><span class="p">({</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">5</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">5</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">5</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;YWIDTH&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">})</span>
        <span class="k">if</span> <span class="n">rotate</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
            <span class="n">parinfo</span><span class="o">.</span><span class="n">append</span><span class="p">({</span><span class="s">&#39;n&#39;</span><span class="p">:</span><span class="mi">6</span><span class="p">,</span><span class="s">&#39;value&#39;</span><span class="p">:</span><span class="n">params</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span><span class="s">&#39;limits&#39;</span><span class="p">:[</span><span class="n">minpars</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span><span class="n">maxpars</span><span class="p">[</span><span class="mi">6</span><span class="p">]],</span><span class="s">&#39;limited&#39;</span><span class="p">:[</span><span class="n">limitedmin</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span><span class="n">limitedmax</span><span class="p">[</span><span class="mi">6</span><span class="p">]],</span><span class="s">&#39;fixed&#39;</span><span class="p">:</span><span class="n">fixed</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span><span class="s">&#39;parname&#39;</span><span class="p">:</span><span class="s">&quot;ROTATION&quot;</span><span class="p">,</span><span class="s">&#39;error&#39;</span><span class="p">:</span><span class="mi">0</span><span class="p">})</span>

    <span class="k">if</span> <span class="n">extra_pars</span><span class="p">:</span>
        <span class="k">for</span> <span class="n">P</span> <span class="ow">in</span> <span class="n">extra_pars</span><span class="p">:</span>
            <span class="n">parinfo</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">P</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">autoderiv</span> <span class="ow">is</span> <span class="bp">False</span><span class="p">:</span>
        <span class="c"># the analytic derivative, while not terribly difficult, is less</span>
        <span class="c"># efficient and useful.  I only bothered putting it here because I was</span>
        <span class="c"># instructed to do so for a class project - please ask if you would</span>
        <span class="c"># like this feature implemented</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&quot;I&#39;m sorry, I haven&#39;t implemented this feature yet.&quot;</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
<span class="c">#        p, cov, infodict, errmsg, success = optimize.leastsq(errorfunction,\</span>
<span class="c">#                params, full_output=1)</span>
        <span class="n">mp</span> <span class="o">=</span> <span class="n">mpfit</span><span class="p">(</span><span class="n">mpfitfun</span><span class="p">(</span><span class="n">data</span><span class="p">,</span><span class="n">err</span><span class="p">),</span><span class="n">parinfo</span><span class="o">=</span><span class="n">parinfo</span><span class="p">,</span><span class="n">quiet</span><span class="o">=</span><span class="n">quiet</span><span class="p">)</span>


    <span class="k">if</span> <span class="n">returnmp</span><span class="p">:</span>
        <span class="n">returns</span> <span class="o">=</span> <span class="p">(</span><span class="n">mp</span><span class="p">)</span>
    <span class="k">elif</span> <span class="n">return_parinfo</span><span class="p">:</span>
        <span class="n">returns</span> <span class="o">=</span> <span class="p">(</span><span class="n">parinfo</span><span class="p">)</span>
    <span class="k">elif</span> <span class="n">return_all</span> <span class="ow">is</span> <span class="bp">False</span><span class="p">:</span>
        <span class="n">returns</span> <span class="o">=</span> <span class="n">mp</span><span class="o">.</span><span class="n">params</span>
    <span class="k">elif</span> <span class="n">return_all</span><span class="p">:</span>
        <span class="n">returns</span> <span class="o">=</span> <span class="n">mp</span><span class="o">.</span><span class="n">params</span><span class="p">,</span><span class="n">mp</span><span class="o">.</span><span class="n">perror</span>
    <span class="k">if</span> <span class="n">returnfitimage</span><span class="p">:</span>
        <span class="n">fitimage</span> <span class="o">=</span> <span class="n">psffunction</span><span class="p">(</span><span class="n">mp</span><span class="o">.</span><span class="n">params</span><span class="p">,</span><span class="n">circle</span><span class="p">,</span><span class="n">rotate</span><span class="p">,</span><span class="n">vheight</span><span class="p">)(</span><span class="o">*</span><span class="n">numpy</span><span class="o">.</span><span class="n">indices</span><span class="p">(</span><span class="n">data</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span>
        <span class="n">returns</span> <span class="o">=</span> <span class="p">(</span><span class="n">returns</span><span class="p">,</span><span class="n">fitimage</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">returns</span></div>
</pre></div>

          </div>
        </div>
      </div>
        </div>
        <div class="sidebar">
          <h3>Table Of Contents</h3>
          <ul>
<li class="toctree-l1"><a class="reference internal" href="../../agpy.html">Adam Ginsburg&#8217;s Python Code (agpy)</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../image_tools.html">Image Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../fft_tools.html">AG_fft_tools Package</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../plfit.html">plfit Package</a></li>
</ul>

          <h3 style="margin-top: 1.5em;">Search</h3>
          <form class="search" action="../../search.html" method="get">
            <input type="text" name="q" />
            <input type="submit" value="Go" />
            <input type="hidden" name="check_keywords" value="yes" />
            <input type="hidden" name="area" value="default" />
          </form>
          <p class="searchtip" style="font-size: 90%">
            Enter search terms or a module, class or function name.
          </p>
        </div>
        <div class="clearer"></div>
      </div>
    </div>

    <div class="footer">
      &copy; Copyright 2011, Adam Ginsburg.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.2pre.
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setAccount', 'UA-6253248-2']);
      _gaq.push(['_trackPageview']);

      (function() {
        var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
        ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
        var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
      })();

    </script>
        
    </div>
  </body>
</html>